function varargout = anncyc (s, ps)

   ## usage:  u = anncyc (s, ps)
   ##
   ## estimate and remove annual cycle, including scaling

   persistent nprint = 0 ;
   global PAR PRJ BATCH

   smin = PAR.scyc.fun ;
   if isempty(smin)
      lp0 = 0 ;
   else
      p0 = PAR.scyc.init ; lp0 = length(p0) ;
   endif

   u = s ;
   N = size(u.x) ;
   j0 = length(N)-1 ;
   K = [1 N(end-j0+1:end)] ;

   if (nargin < 2)
      nprint++ ;
      ps = nan(PAR.numhar+lp0, prod(K)) ;
   else
      ps = reshape(ps, size(ps,1), prod(K(2:end))) ;
   endif 
   u.x = reshape(u.x, prod(N(1:end-j0)), prod(K)) ;
   if !isfield(s, "cal"), s.cal = "gregorian" ; endif 
   t = date2toy(u.id, s.cal) ;
   if (PAR.numhar == 3)
      x_h = [ones(rows(t),1), cos(t), sin(t)] ;
   elseif (PAR.numhar == 5)
      x_h = [ones(rows(t),1) cos(t) sin(t) cos(2*t) sin(2*t)] ;
   endif 

   if (nargin < 2)
      CJ = mat2cell(u.x, N(1), ones(1, N(2))) ;
      ps = parcellfun_(@(x) _regress(x, x_h), CJ, PAR.par_opt{:}) ;
      ps = cell2mat(ps) ;

      if !isempty(smin) ## Taiga workaround
	 u.a = u.x - x_h * ps ;
	 u.a(u.a == -inf) = nan ; # deal with ispos
	 CJ = mat2cell(abs(u.a), N(1), ones(1, N(2))) ;
	 ps2 = parcellfun_(@(x) fminunc(@(p) sumsq(x-smin(p,t)), p0), CJ, PAR.par_opt{:}) ;
	 ps2 = cell2mat(ps2) ;
	 ps = [ps ; ps2] ;

	 if PAR.dbg
	    if BATCH
	       graphics_toolkit gnuplot ; figure("visible", "off") ;
	    endif
	    subplot(5,1,1) ; plot(pltid(u.id), [u.x x_h * ps(1:PAR.numhar,:)]) ;	    
	    set(gca, "colororder", jet(columns(ps2))) ;
	    subplot(5,1,2) ; plot(pltid(u.id), u.a) ;
	    set(gca, "colororder", jet(columns(ps2))) ;
	    subplot(5,1,3) ; plot(pltid(u.id), [abs(u.a) smin(ps2,t)]) ;
	    set(gca, "colororder", jet(columns(ps2))) ;
	    subplot(5,1,4) ; plot(pltid(u.id), [abs(u.a)./smin(ps2,t)]) ;
	    set(gca, "colororder", jet(columns(ps2))) ;
	    subplot(5,1,5) ; plot(pltid(u.id), [u.a ./ smin(ps2,t)]) ;
	    set(gca, "colororder", jet(columns(ps2))) ;
	    printf("anncyc: --> %s\n", pfile = fullfile("data", PRJ, PAR.pdd, ["anc_" num2str(nprint) ".png"])) ;
	    print(pfile) ; clf ;

	 endif
      endif

      if (nargout < 2)
	 varargout{1} = reshape(ps, [rows(ps) K(2:end)]) ;
	 return ;
      endif
   endif 

   u.a = u.x - x_h * ps(1:PAR.numhar,:) ;

   if !isempty(smin)
      ps2 = ps(PAR.numhar+1:end,:) ;
      u.a = u.a ./ smin(ps2,t) ;
   endif

   u.x = reshape(u.x, N) ;
   u.a = reshape(u.a, N) ;

   if (nargin > 1)
      varargout{1} = u ;
   else
      varargout{1} = reshape(ps, [rows(ps) K(2:end)]) ;
      varargout{2} = u ;
   endif 
      
endfunction


function res = _regress (y, x)

   ## usage:  res = regresswr (y, x)
   ##
   ## wrapper for regress

   global PAR ;

   I = !isinf(y) ;
   if (sum(!isnan(y)) < PAR.Ccases)
      res = nan(5,1) ;
   else
      res = regress(y(I), x(I,:)) ;
   endif 

endfunction
